Diversity patterns of the Flora of Greece with applications in R
Όλα τα απαραίτητα αρχεία βρίσκονται στο github.
1 Πριν ξεκινήσουμε
Βεβαιώστε ότι το όνομα του project σας, όπως και το file path ΔΕΝ περιέχουν:
1. Ελληνικούς χαρακτήρες
2. κενά (αντικαταστήστε τα κενά με _)
Επί παραδείγματι το file path σας θα μπορούσε να είναι: C:/Environmental_Data_Analysis/Geographical_Data.
Όλα τα απαιτούμενα δεδομένα και αρχεία είναι ελεύθερα διαθέσιμα στο διαδίκτυο.
2 Φόρτωση βιβλιοθηκών
Ας φορτώσουμε όλες τις βιβλιοθήκες τις οποίες θα χρησιμοποιήσουμε σήμερα1.
## ===========================================================================##
## Load them
## ===========================================================================##
library(raster)
library(sf)
library(tidyverse)
library(ConR)
library(red)
library(rasterVis)
library(ggspatial)
library(picante)
library(caper)
library(phyloregion)
library(ggsflabel)
library(cowplot)
## ===========================================================================##3 Διανυσματικά δεδομένα
3.1 Σύνορα χωρών
Θα χρησιμοποιήσουμε τα διανυσματικά δεδομένα της Ιταλίας, όπως έχουμε κάνει και σε προηγούμενα εργαστήρια.
Θυμηθείτε ότι θα χρειαστεί να αλλάξετε το file path και να χρησιμοποιήσετε τα δεδομένα από το Github.
## ===========================================================================##
## Load the spatial data for Italian provinces
## ===========================================================================##
Italy <- read_sf("Italin provinces.shp") %>% st_transform(4326) %>% as_Spatial()
## ===========================================================================##
## ===========================================================================##
## Load the spatial data for Italian counties
## ===========================================================================##
Italian_counties <- read_sf("reg2011_g.shp") %>% mutate(NOME_REG = str_replace_all(NOME_REG,
c(`VALLE D'AOSTA/VALLÉE D'AOSTE\r\nVALLE D'AOSTA/VALLÉE D'AOSTE` = "AOSTA", `TRENTINO-ALTO ADIGE/SUDTIROL` = "TRENTINO",
`FRIULI VENEZIA GIULIA` = "FRIULI", `EMILIA-ROMAGNA` = "EMILIA"))) %>% st_transform(4326)
## ===========================================================================##
## ===========================================================================##
## Keep only the names in Trentino
## ===========================================================================##
Italy_sf_tr <- Italy %>% st_as_sf() %>% st_join(Italian_counties) %>% filter(str_detect(NOME_REG,
"TREN"))
## ===========================================================================##4 Δεδομένα θέσης
4.1 Διαδικτυακή βάση δεδομένων GBIF
Θα χρησιμοποιήσουμε τα δεδομένα θέσης για όλα τα φυτικά είδη τα οποία απαντώνται στην Ιταλία, τα οποία περιέχονται στην διαδικτυακή βάση δεδομένων GBIF και από τα οποία έχουμε αφαιρέσει τις όποιες προβληματικές καταγραφές. Θα χρειαστεί να καθαρίσουμε τα ονόματα των taxa και στη συνέχεια να κρατήσουμε μόνο τα έγκυρα ονόματα.
## ============================================================================##
## Load the biotic cleaned data
## ============================================================================##
data_it_cleaned <- readRDS("Italian data cleaned.rds")
## ============================================================================##
## ============================================================================##
## First load the excel file with the correct species names
## ============================================================================##
nms <- readxl::read_excel("italian_names.xlsx")
## ============================================================================##
## ============================================================================##
## Keep only one instance per taxon from our occurrence data
## ============================================================================##
nms_flags <- data_it_cleaned %>% distinct(scientificName) %>% mutate(Taxon = nms$scientificName,
action = nms$action) %>% as_tibble()
## ============================================================================##
## ============================================================================##
## Then create the coords for the correct taxa
## ============================================================================##
coords <- data_it_cleaned %>% full_join(nms_flags) %>% as_tibble() %>% filter(!str_detect(action,
"DEL")) %>% dplyr::select(decimalLongitude, decimalLatitude, Taxon)
clean_dts <- data_it_cleaned %>% full_join(nms_flags) %>% as_tibble() %>% filter(!str_detect(action,
"DEL")) %>% distinct(Taxon, .keep_all = T)
## ============================================================================##5 Αξιολόγηση κατά IUCN
Η αξιολόγηση του κινδύνου εξαφάνισης ενός είδους βασίζεται στην τυποποιημένη διαδικασία που έχει αναπτυχθεί από την Διεθνή Ένωση για την Προστασία της Φύσης (IUCN) και αποτελεί την πλέον αντικειμενική και αναλυτική προσέγγιση για την ανάδειξη αξόνων προτεραιότητας για την διατήρηση των ειδών και οικοτόπων. Ο Ερυθρός Κατάλογος των Απειλούμενων Ειδών ή IUCN Red List of Threatened Species παρέχει πληροφορίες σχετικά με την ταξινόμηση, την κατανομή, καθώς και το καθεστώς διατήρησης των φυτικών και ζωικών ειδών του πλανήτη μας, βάσει των κριτηρίων και των κατηγοριών της IUCN (IUCN Red List Categories and Criteria). Ο κύριος στόχος της διαδικασίας αυτής είναι να προσδιορίσει τον ενδεχόμενο κίνδυνο εξαφάνισης ενός είδους και κατά συνέπεια, να αναδείξει τα είδη εκείνα τα οποία χρήζουν άμεσης προτεραιότητας για την διατήρηση τους (Stévart et al., 2019).
Η δημιουργία ενός Ερυθρού Καταλόγου απαιτεί αξιόπιστα δεδομένα και προσεκτική εφαρμογή όλων των κριτηρίων της IUCN – μια διαδικασία η οποία είναι ιδιαίτερα χρονοβόρα. Το γεγονός αυτό αντανακλάται από το ότι υπάρχουν σημαντικά κενά στα έως τώρα αξιολογηθέντα είδη από την IUCN, καθώς μόλις 28.114 φυτικά είδη έχουν αξιολογηθεί μέχρι σήμερα. Το κενό αυτό είναι απόρροια της εξαιρετικά μεγάλης ποικιλότητας των φυτικών ειδών. Συνεπώς, παρότι η IUCN είναι κοντά στον στόχο που είχε θέσει, να έχει αξιολογήσει δηλαδή 38.500 φυτικά είδη έως το 2020 (Bachman et al., 2018), εντούτοις, ο αριθμός αυτός αποτελεί μόλις το 10% της φυτικής ποικιλότητας του πλανήτη και του δεύτερου στόχου της Παγκόσμιας Στρατηγικής για την Διατήρηση της Φυτικής ποικιλότητας (i.e., Global Strategy for Plant Conservation - περισσότερες πληροφορίες εδώ). Το γεγονός αυτό ώθησε την επιστημονική κοινότητα να προσπαθήσει να βρει τρόπους ώστε να επιταχύνει την διαδικασία αυτή.
6 Εκτίμηση κινδύνου εξαφάνισης - Εισαγωγικές έννοιες
Όπως είπαμε, μέχρι στιγμής έχει αξιολογηθεί ένα πολύ μικρό ποσοστό της έως τώρα γνωστής βιοποικιλότητας:
Η πλειονότητα των χωρών που αντιμετωπίζουν τα μεγαλύτερα προβλήματα όσον αφορά την κατάσταση διατήρησης των φυτικών taxa τα οποία φιλοξενούν εντοπίζονται στις λιγότερο αναπτυγμένες χώρες του πλανήτη:
Σε ευρωπαϊκό επίπεδο δε, γνωρίζουμε ότι η κατάσταση έχει προς το παρόν διαμορφωθεί ως εξής:
Όπου μπορούμε να διαπιστώσουμε ότι οι χώρες της Μεσογειακής λεκάνης και ιδιαίτερα η Ελλάδα, θεωρούνται θερμά σημεία κινδύνου εξαφάνισης:
Με την κυριότερη απειλή να είναι η εντατικοποίηση της βόσκησης
Όμως οι απειλές αυτές διαφέρουν τόσο μεταξύ των θερμών και των ψυχρών σημείων βιοποικιλότητας:
Όσο και στον χρόνο (για τα τελευταία 250 έτη τουλάχιστον):
Σύμφωνα με την IUCN, υπάρχουν 5 κατηγορίες κινδύνου εξαφάνισης2:
Ας τα δούμε πιο αναλυτικά:
Όσον αφορά δε την Ελλάδα:
Ενώ όσον αφορά την Κρήτη, το θερμότερο σημείο ενδημικής φυτικής ποικιλότητας στη Μεσόγειο:
Η κλιματική αλλαγή αναμένεται να έχει μεγάλες επιπτώσεις όσον αφορά τα θερμά σημεία ενδημικής φυτικής ποικιλότητας της νήσου (Kougioumoutzis et al., 2020a, 2020b)
7 Εκτίμηση κινδύνου εξαφάνισης - Στην πράξη
7.1 Φόρτωση των δεδομένων θέσης
Τώρα θα φορτώσουμε το set δεδομένων μας.
Στην περίπτωση που δεν έχετε κατεβάσει ήδη τα απαραίτητα αρχεία, κάντε το τώρα:
## ==================================================================================##
## Load the data
## ==================================================================================##
cretan_taxa <- readxl::read_excel("./PVA_Crete.xlsx") %>% dplyr::rename(decimalLongitude = lon,
decimalLatitude = lat)
cretan_taxa$Taxon <- factor(cretan_taxa$Taxon)
levels(cretan_taxa$Taxon)## [1] "Allium dilatatum"
## [2] "Campanula jacquinii subsp. jacquinii"
## [3] "Tulipa doerfleri"
7.2 Μορφοποιήση των δεδομένων
Το πρώτο πράγμα το οποίο θα χρειαστεί να κάνουμε, είναι να μορφοποιήσουμε με τον κατάλληλο τρόπο τα δεδομένα μας, ώστε να μπορέσουμε να τρέξουμε τις εντολές από την βιβλιοθήκη red. Θα χρειαστεί να φτιάξουμε τρια αντικείμενα, ένα για κάθε είδος το οποίο μας ενδιαφέρει και για το οποίο έχουμε δεδομένα. Σήμερα θα ασχοληθούμε τρια από τα τοπικά νησιωτικά ενδημικά της Κρήτης, το Allium dilalatum, την Campanula jacquinii subsp. jacquinii και την Tulipa doerfleri.
Το είδος Allium dilalatum - Photo from www.greekflora.gr
Το είδος Campanula jacquinii subsp. jacquinii - Photo from www.west-crete.com
Το είδος Tulipa doerfleri- Photo from idr
## ==================================================================================##
## Create 3 objects
## ==================================================================================##
all_dila <- cretan_taxa %>% filter(Taxon == "Allium dilatatum")
camp_jac <- cretan_taxa %>% filter(Taxon == "Campanula jacquinii subsp. jacquinii")
tul_doer <- cretan_taxa %>% filter(Taxon == "Tulipa doerfleri")
## ==================================================================================##7.3 Μετατροπή σε μήτρα (matrix)
Πρέπει να μετατρέψουμε τα ανωτέρω αντικείμενα σε μήτρα προκειμένου να είναι συμβατά με την βιβλιοθήκη red:
## ==================================================================================##
## Convert to matrices
## ==================================================================================##
all_dila_mat <- all_dila %>% dplyr::select(decimalLongitude:decimalLatitude) %>%
as.matrix()
camp_jac_mat <- camp_jac %>% dplyr::select(decimalLongitude:decimalLatitude) %>%
as.matrix()
tul_doer_mat <- tul_doer %>% dplyr::select(decimalLongitude:decimalLatitude) %>%
as.matrix()
## ==================================================================================##7.4 Εκτίμηση περιοχής κατάληψης (Area of Occupancy)
Τώρα είμαστε σε θέση να υπολογίσουμε την παρούσα περιοχή κατάληψης των ειδών που μας ενδιαφέρουν:
## ==================================================================================##
## Area of occupancy
## ==================================================================================##
aoo(all_dila_mat)
aoo(camp_jac_mat)
aoo(tul_doer_mat)
## ==================================================================================##7.5 Εκτίμηση περιοχής εμφάνισης (Extent of Occurrence)
Τώρα είμαστε σε θέση να υπολογίσουμε την παρούσα περιοχή εμφάνισης των ειδών που μας ενδιαφέρουν:
## ==================================================================================##
## Extent of Occurrence
## ==================================================================================##
eoo(all_dila_mat)
eoo(camp_jac_mat)
eoo(tul_doer_mat)
## ==================================================================================##7.6 Εκτίμηση περιοχής κατάληψης (Area of Occupancy) στο μέλλον
Προκειμένου να υπολογίσουμε την παρούσα περιοχή κατάληψης των ειδών που μας ενδιαφέρουν, θα χρειαστεί να φορτώσουμε την πρόβλεψη για την κατάσταση τους στο μέλλον και δη για το έτος 2070:
## ==================================================================================##
## Load the rds files
## ==================================================================================##
allium_future <- readRDS("./RDS/Allium future.rds")
campanula_future <- readRDS("./RDS/Campanula future.rds")
tulipa_future <- readRDS("./RDS/Tulipa future.rds")
## ==================================================================================##
## ==================================================================================##
## Area of occupancy
## ==================================================================================##
aoo(allium_future)
aoo(campanula_future)
aoo(tulipa_future)
## ==================================================================================##7.7 Εκτίμηση περιοχής εμφάνισης (Extent of Occurrence) στο μέλλον
Ας υπολογίσουμε την παρούσα περιοχή εμφάνισης των ειδών που μας ενδιαφέρουν:
## ==================================================================================##
## Extent of Occurrence
## ==================================================================================##
eoo(allium_future)
eoo(campanula_future)
eoo(tulipa_future)
## ==================================================================================##Τι παρατηρείτε; Υπάρχει διαφορά σε σχέση με την παρούσα κατάσταση; Πώς θα χαρακτηρίζατε τα είδη αυτά σύμφωνα με τις κατηγορίες κινδύνου εξαφάνισης της IUCN;
8 Ταχεία και μαζική αξιολόγηση ειδών κατά IUCN – The PACA approach
Η αξιολόγηση του καθεστώτος κινδύνου εξαφάνισης ενός είδους βάσει των κατευθυντήριων γραμμών της Διεθνούς Ενώσεως για την Διατήρηση της Φύσης IUCN αποτελούν τον πλέον αντικειμενικό και αναλυτικό τρόπο προσέγγισης για τον καθορισμό των στόχων προτεραιότητας διατήρησης. Ο ιστότοπος IUCN Red List of Threatened Species παρέχει πλειάδα πληροφοριών σχετικά με την ταξινομική, την κατανομή και το καθεστώς διατήρησης αρκετών φυτικών και ζωικών ειδών, βάσει των Κατηγοριών και των Κριτηρίων της IUCN. Ο κύριος στόχος αυτής της διαδικασίας είναι να προσδιοριστεί, χρησιμοποιώντας μια αναλυτικότατη και αντικειμενική μέθοδο, ο κίνδυνος εξαφάνισης ενός είδους και, συνεπώς, να προσδιοριστεί ποια είδη χρήζουν άμεσης προτεραιότητας διατήρησης. Ο κατάλογος Ερυθρών Δεδομένων της IUCN θεωρείται ως ο “χρυσός κανόνας” όσον αφορά τον εντοπισμό των ειδών εκείνων για τα οποία απαιτείται να δοθεί ιδιαίτερη προσοχή κατά τον σχεδιασμό και την υλοποίηση έργων περιβαλλοντικών επιπτώσεων (Stévart et al., 2019).
Η αξιολόγηση συνεπώς του καθεστώτος κινδύνου εξαφάνισης βάσει των Κριτηρίων της IUCN απαιτεί τόσο αξιόπιστα δεδομένα, όσο και προσεκτικότατη εφαρμογή και ενδελεχή γνώση των Κριτηρίων της IUCN και για αυτόν τον λόγο είναι μια ιδιαίτερα χρονοβόρα διαδικασία. Αυτό υποδηλώνεται και από το γεγονός ότι ένα μικρό ποσοστό των φυτικών ειδών (~10%) έχει έως τώρα αξιολογηθεί σύμφωνα με τα Κριτήρια της IUCN (Bachman et al., 2018) και από ότι φαίνεται ο Aichi Target 2 (Αξιολόγηση του καθεστώτος κινδύνου εξαφάνισης για την πλειονότητα των φυτικών ειδών - Global Strategy for Plant Conservation και εδώ), δεν πρόκειται να επιτευχθεί σύντομα. Για αυτόν τον λόγο, η επιστημονική κοινότητα καλέιται να επισπεύσει - στο μέτρο του δυνατού - την διαδικασία αξιολόγησης των φυτικών ειδών σε παγκόσμιο επίπεδο.
Πρόσφατα, προτάθηκε μια νέα προσέγγιση που βασίζεται στα κριτήρια της IUCN και παρέχει την δυνατότητα ταχείας, μαζικής και αξιόπιστης αξιολόγησης των ειδών μιας περιοχής/χώρας (Dauby et al., 2017; Stévart et al., 2019). Η προσέγγιση αυτή ονομάζεται Preliminary Automated Conservation Assessments (PACA) και αυτή θα χρησιμοποιήσουμε σήμερα.
Υπάρχει πιθανότητα αναλόγως των δυνατοτήτων του Η/Υ σας, η εντολή αυτή να διαρκέσει μεγάλο χρονικό διάστημα. Μπορείτε να εισαγάγετε το αποτέλεσμα αυτής με την εντολή:
italy_iucn_eval <- readRDS('Italian IUCN.rds').
##============================================================================##
## First, let's select only the required columns
##============================================================================##
data_it_sel <- coords %>%
dplyr::select(decimalLatitude,
decimalLongitude,
Taxon)
##============================================================================##
##============================================================================##
## Then, let's evaluate our taxa
##============================================================================##
italy_iucn_eval <- IUCN.eval(data_it_sel,
## Here we use the spatial object for the country
## we want
country_map = Italy,
## We can run the function in parallel to save
## time
parallel = T,
## Change the number of cores accordingly
NbeCores = 23 ) %>%
## Here we create a new column (EOO_tr) and we substitute the NA values with
## the values from the AOO column (Why?)
mutate(EOO_tr = ifelse(is.na(EOO),
AOO,
EOO)) %>%
as_tibble() %>%
dplyr::select(taxa,
EOO_tr,
everything())
##============================================================================##
## Let's see the result
##============================================================================##
italy_iucn_evalΜε αυτόν τον τρόπο λοιπόν, αξιολογήσαμε το καθεστώς κινδύνου εξαφάνισης όλων των ειδών που είχαμε κατεβάσει από την GBIF για την Ιταλία. Όμως, παρότι είχαμε προχωρήσει σε ορισμένα βήματα καθαρισμού των δεδομένων, απαιτούνται πολύ περισσότερα στάδια3 προκειμένου να αξιολογηθεί ένα είδος σύμφωνα με όλα τα κριτήρια της IUCN.
Μπορούμε επίσης να υπολογίσουμε και να αποθηκεύσουμε στον σκληρό μας δίσκο τα πολύγωνα της περιοχής εμφάνισης (Extent of Occurrence) για κάθε είδος, βασιζόμενοι σε δύο μεθόδους προσδιορισμού του πολυγώνου αυτού: την convex- και την alpha-hull μέθοδο (Εάν θέλετε, μπορείτε να διαβάσετε αυτήν την εργασία, εάν ενδιαφέρεστε να μάθετε κάτι περισσότερο σχετικά με τις μεθόδους αυτές).
Μην τρέξετε αυτή την εντολή τώρα.
##============================================================================##
## First we run the convex-hull method
##============================================================================##
eooshp_convexhull <- EOO.computing(data_gr_sel,
exclude.area = T,
country_map = Greece,
export_shp = T,
write_shp = T,
Name_Sp = 'scientificName',
parallel = T,
NbeCores = 4, ## Change them
show_progress = T,
method.less.than3 = 'arbitrary') ## What does this option mean?
convexhull <- eooshp_convexhull$spatial.polygon_1
##============================================================================##
##============================================================================##
## Then the alpha-hull method
##============================================================================##
eooshp_ahull <- EOO.computing(data_gr_sel,
exclude.area = T,
country_map = Greece,
export_shp = T,
write_shp = T,
Name_Sp = 'scientificName',
parallel = T,
NbeCores = 4, ## Change them
show_progress = T,
method.less.than3 = 'arbitrary',
method.range = 'alpha.hull')
ahull <- eooshp_ahull$spatial.polygon_1
##============================================================================##8.1 Χωρικά πρότυπα του κινδύνου εξαφάνισης
Ας δούμε τώρα ποια περιοχή της Ιταλίας εμφανίζει το υψηλότερο ποσοστό κινδυνευόντων φυτικών ειδών:
##============================================================================##
## First create a common column and then join the occurrences with the
## evaluations
##============================================================================##
italy_iucn_eval <- italy_iucn_eval %>% rename(Taxon = taxa)
data_joined <- full_join(coords, italy_iucn_eval)
##============================================================================##
##============================================================================##
## Then select only the endangered taxa (CR, EN, VU)
##============================================================================##
endangered <- data_joined %>% filter(!str_detect(Category_code, 'LC'))
##============================================================================##
##============================================================================##
## Convert Italy to an sf spatial object
##============================================================================##
Italy_sf <- Italy %>% st_as_sf()
##============================================================================##
##============================================================================##
## Convert the endangered species to a sf spatial object
##============================================================================##
data_it_sf <- st_as_sf(endangered,
coords = c("decimalLongitude",
"decimalLatitude"),
crs = 4326)
##============================================================================##
##============================================================================##
## Make sure that all these occurrences lie within Greece
##============================================================================##
point_endangered <- data_it_sf[Italy_sf,]
##============================================================================##
##============================================================================##
## Intersect Greece with these occurrences
##============================================================================##
endangered_intersection <- Italy_sf %>%
st_join(point_endangered) %>%
group_by(NOME_PRO) %>% ## Here we group the data based on the county names
## First we create a variable containing the sum of the species in each county
summarize(Species_Number = sum(NROW(unique(Taxon))),
## Then we just sum the number of occurrences in each county
Occurrence_Number = n(),
## Finally, we calculate the proportion of the endangered species in
## each county
Risk = sum(NROW(unique(Taxon)))/sum(NROW(unique(data_joined$Taxon)))*100
)
##============================================================================##
##============================================================================##
## Colour function
##============================================================================##
col2alpha <- function(col, alpha) {
col_rgb <- col2rgb(col)/255
rgb(col_rgb[1], col_rgb[2], col_rgb[3],
alpha = alpha)
}
##============================================================================####============================================================================##
## Create the main map
##============================================================================##
main_map <- ggplot(Italy_sf) +
geom_sf(fill = "antiquewhite1") +
geom_sf() +
geom_sf(data = Italy_sf,
fill = NA) +
geom_sf(aes(fill = Risk),
data = endangered_intersection,
inherit.aes = FALSE) +
theme(axis.title = element_blank(),
legend.position = "bottom") +
coord_sf(expand = FALSE) +
scale_fill_viridis_c(option = "viridis",
trans = "sqrt",
limits = c(0, 12.4),
breaks = c(0.4, 3, 12.3),
labels = scales::label_percent(scale = 1),
name = expression(bold('Extinction risk'))) +
geom_rect(
xmin = st_bbox(Italy_sf_tr)[[1]],
ymin = st_bbox(Italy_sf_tr)[[2]],
xmax = st_bbox(Italy_sf_tr)[[3]],
ymax = st_bbox(Italy_sf_tr)[[4]],
fill = NA,
colour = "red",
size = 0.6
) +
annotation_scale(location = "bl",
width_hint = 0.25) +
annotation_north_arrow(location = "bl",
which_north = "true",
pad_x = unit(0.75, "in"),
pad_y = unit(0.5, "in"),
style = north_arrow_fancy_orienteering) +
geom_sf(data = Italian_counties,
fill = NA,
lwd = 1,
color = 'black') +
geom_sf_label_repel(aes(label = NOME_REG),
data = Italian_counties) +
theme(panel.grid.major = element_line(color = gray(.5),
linetype = "dashed",
size = 0.5),
panel.background = element_rect(fill = col2alpha('steelblue', 0.2))) +
labs(title = 'Extinction risk in Italian counties',
subtitle = 'Source: GBIF')
##============================================================================##
##============================================================================##
## Create the inset map
##============================================================================##
side_map <- ggplot(Italy_sf) +
geom_sf(fill = "antiquewhite1") +
geom_sf() +
geom_sf(data = Italy_sf,
fill = NA) +
geom_sf(aes(fill = Risk),
data = endangered_intersection,
inherit.aes = FALSE) +
theme(axis.title = element_blank(),
legend.position = "bottom") +
coord_sf(expand = FALSE) +
scale_fill_viridis_c(option = "viridis",
trans = "sqrt",
limits = c(0, 12.4),
breaks = c(0.4, 3, 12.3),
labels = scales::label_percent(scale = 1),
name = 'Extinction risk') +
geom_rect(
xmin = st_bbox(Italy_sf_tr)[[1]],
ymin = st_bbox(Italy_sf_tr)[[2]],
xmax = st_bbox(Italy_sf_tr)[[3]],
ymax = st_bbox(Italy_sf_tr)[[4]],
fill = NA,
colour = "red",
size = 0.6
) +
geom_sf(data = Italian_counties,
fill = NA,
lwd = 1,
color = 'black') +
geom_sf_label_repel(aes(label = NOME_PRO),
data = Italy_sf_tr) +
coord_sf(xlim = c(st_bbox(Italy_sf_tr)[[1]],
st_bbox(Italy_sf_tr)[[3]]),
ylim = c(st_bbox(Italy_sf_tr)[[2]],
st_bbox(Italy_sf_tr)[[4]]),
expand = FALSE) +
theme_void() +
theme(legend.position = "none")
##============================================================================##
##============================================================================##
## And the final map
##============================================================================##
main_map %>%
ggdraw() +
draw_plot(side_map,
# The distance along a (0,1) x-axis to draw the left edge of the plot
x = 0.625,
# The distance along a (0,1) y-axis to draw the bottom edge of the plot
y = 0.65,
# The width and height of the plot expressed as proportion of the entire
# ggdraw object
width = 0.35,
height = 0.35)Ποια περιοχή της Ιταλίας εμφανίζει το υψηλότερο ποσοστό κινδυνευόντων ειδών;
## ============================================================================##
## The solution
## ============================================================================##
endangered_intersection %>% arrange(desc(Risk))| NOME_PRO | Species_Number | Occurrence_Number | Risk | geometry |
|---|---|---|---|---|
| BOLZANO/BOZEN | 783 | 2108 | 12.3990499 | MULTIPOLYGON (((12.20557 47… |
| COMO | 467 | 2585 | 7.3950911 | MULTIPOLYGON (((9.279479 46… |
| VALLE D’AOSTA | 400 | 1020 | 6.3341251 | MULTIPOLYGON (((7.58858 45…. |
| VERBANO-CUSIO-OSSOLA | 392 | 1409 | 6.2074426 | MULTIPOLYGON (((8.449767 46… |
| SONDRIO | 347 | 1118 | 5.4948535 | MULTIPOLYGON (((10.2493 46…. |
| L’AQUILA | 345 | 695 | 5.4631829 | MULTIPOLYGON (((13.40447 42… |
| TORINO | 302 | 467 | 4.7822644 | MULTIPOLYGON (((7.859052 45… |
| TRENTO | 279 | 413 | 4.4180523 | MULTIPOLYGON (((11.82411 46… |
| CUNEO | 267 | 495 | 4.2280285 | MULTIPOLYGON (((7.990904 44… |
| VARESE | 266 | 884 | 4.2121932 | MULTIPOLYGON (((8.78101 46…. |
| PALERMO | 264 | 403 | 4.1805226 | MULTIPOLYGON (((13.07012 38… |
| GROSSETO | 214 | 639 | 3.3887569 | MULTIPOLYGON (((11.13574 42… |
| BERGAMO | 199 | 301 | 3.1512272 | MULTIPOLYGON (((10.10323 46… |
| BELLUNO | 184 | 260 | 2.9136975 | MULTIPOLYGON (((12.50593 46… |
| FOGGIA | 169 | 392 | 2.6761679 | MULTIPOLYGON (((16.06218 41… |
| UDINE | 166 | 250 | 2.6286619 | MULTIPOLYGON (((12.76338 46… |
| NUORO | 157 | 283 | 2.4861441 | MULTIPOLYGON (((9.583129 40… |
| BRESCIA | 150 | 188 | 2.3752969 | MULTIPOLYGON (((10.49807 46… |
| SIRACUSA | 150 | 293 | 2.3752969 | MULTIPOLYGON (((15.11665 37… |
| ALESSANDRIA | 138 | 222 | 2.1852732 | MULTIPOLYGON (((8.404617 45… |
| CATANIA | 138 | 270 | 2.1852732 | MULTIPOLYGON (((15.09013 37… |
| OGLIASTRA | 136 | 198 | 2.1536025 | MULTIPOLYGON (((9.628694 40… |
| TRAPANI | 128 | 175 | 2.0269200 | MULTIPOLYGON (((12.44316 37… |
| REGGIO DI CALABRIA | 112 | 136 | 1.7735550 | MULTIPOLYGON (((15.72781 37… |
| MESSINA | 109 | 135 | 1.7260491 | MULTIPOLYGON (((14.96671 38… |
| SIENA | 104 | 154 | 1.6468725 | MULTIPOLYGON (((11.42858 43… |
| ROMA | 101 | 160 | 1.5993666 | MULTIPOLYGON (((12.62241 41… |
| FROSINONE | 97 | 178 | 1.5360253 | MULTIPOLYGON (((13.33753 41… |
| BOLOGNA | 95 | 144 | 1.5043547 | MULTIPOLYGON (((11.29036 44… |
| TREVISO | 91 | 106 | 1.4410135 | MULTIPOLYGON (((12.34077 46… |
| GENOVA | 90 | 125 | 1.4251781 | MULTIPOLYGON (((8.686318 44… |
| VERONA | 84 | 110 | 1.3301663 | MULTIPOLYGON (((10.88381 45… |
| VERCELLI | 83 | 158 | 1.3143310 | MULTIPOLYGON (((8.204472 45… |
| PERUGIA | 82 | 149 | 1.2984956 | MULTIPOLYGON (((12.36166 43… |
| SASSARI | 80 | 133 | 1.2668250 | MULTIPOLYGON (((8.203562 40… |
| LIVORNO | 78 | 105 | 1.2351544 | MULTIPOLYGON (((10.54294 42… |
| RIETI | 76 | 102 | 1.2034838 | MULTIPOLYGON (((13.28007 42… |
| IMPERIA | 75 | 114 | 1.1876485 | MULTIPOLYGON (((7.760733 43… |
| OLBIA-TEMPIO | 75 | 112 | 1.1876485 | MULTIPOLYGON (((9.515061 40… |
| TRIESTE | 75 | 101 | 1.1876485 | MULTIPOLYGON (((13.69162 45… |
| CAGLIARI | 70 | 81 | 1.1084719 | MULTIPOLYGON (((9.309161 39… |
| MILANO | 67 | 99 | 1.0609660 | MULTIPOLYGON (((8.93997 45…. |
| PESCARA | 65 | 78 | 1.0292953 | MULTIPOLYGON (((14.07289 42… |
| CHIETI | 64 | 89 | 1.0134600 | MULTIPOLYGON (((14.30305 42… |
| FIRENZE | 64 | 77 | 1.0134600 | MULTIPOLYGON (((11.44765 44… |
| CALTANISSETTA | 62 | 98 | 0.9817894 | MULTIPOLYGON (((13.81375 37… |
| CARBONIA-IGLESIAS | 62 | 104 | 0.9817894 | MULTIPOLYGON (((8.361779 39… |
| SAVONA | 62 | 79 | 0.9817894 | MULTIPOLYGON (((8.192667 44… |
| LA SPEZIA | 61 | 82 | 0.9659541 | MULTIPOLYGON (((9.567304 44… |
| LECCO | 60 | 97 | 0.9501188 | MULTIPOLYGON (((9.415008 46… |
| TARANTO | 60 | 122 | 0.9501188 | MULTIPOLYGON (((17.29291 40… |
| POTENZA | 59 | 85 | 0.9342835 | MULTIPOLYGON (((15.72353 39… |
| LECCE | 57 | 99 | 0.9026128 | MULTIPOLYGON (((18.27608 39… |
| LUCCA | 55 | 74 | 0.8709422 | MULTIPOLYGON (((10.16621 43… |
| AGRIGENTO | 48 | 89 | 0.7600950 | MULTIPOLYGON (((13.02734 37… |
| ORISTANO | 48 | 68 | 0.7600950 | MULTIPOLYGON (((8.497556 40… |
| TERAMO | 48 | 84 | 0.7600950 | MULTIPOLYGON (((13.92029 42… |
| BRINDISI | 47 | 59 | 0.7442597 | MULTIPOLYGON (((17.96702 40… |
| MATERA | 47 | 65 | 0.7442597 | MULTIPOLYGON (((16.2523 40…. |
| PADOVA | 47 | 54 | 0.7442597 | MULTIPOLYGON (((11.83986 45… |
| VICENZA | 47 | 87 | 0.7442597 | MULTIPOLYGON (((11.53894 46… |
| PARMA | 43 | 53 | 0.6809184 | MULTIPOLYGON (((10.1873 45…. |
| COSENZA | 42 | 50 | 0.6650831 | MULTIPOLYGON (((15.81585 39… |
| RAGUSA | 42 | 73 | 0.6650831 | MULTIPOLYGON (((14.49321 36… |
| VENEZIA | 42 | 43 | 0.6650831 | MULTIPOLYGON (((12.79997 45… |
| AVELLINO | 40 | 49 | 0.6334125 | MULTIPOLYGON (((15.19857 41… |
| BARI | 40 | 57 | 0.6334125 | MULTIPOLYGON (((16.65049 41… |
| ISERNIA | 39 | 83 | 0.6175772 | MULTIPOLYGON (((14.30542 41… |
| ASTI | 38 | 40 | 0.6017419 | MULTIPOLYGON (((8.046812 45… |
| SALERNO | 38 | 42 | 0.6017419 | MULTIPOLYGON (((15.32929 40… |
| FORLI’-CESENA | 37 | 44 | 0.5859066 | MULTIPOLYGON (((12.07366 44… |
| MACERATA | 37 | 54 | 0.5859066 | MULTIPOLYGON (((13.72925 43… |
| MEDIO CAMPIDANO | 37 | 44 | 0.5859066 | MULTIPOLYGON (((8.950374 39… |
| NAPOLI | 35 | 40 | 0.5542359 | MULTIPOLYGON (((14.03959 40… |
| RAVENNA | 35 | 47 | 0.5542359 | MULTIPOLYGON (((12.28071 44… |
| PAVIA | 34 | 40 | 0.5384006 | MULTIPOLYGON (((8.847651 45… |
| VITERBO | 34 | 56 | 0.5384006 | MULTIPOLYGON (((11.86313 42… |
| PISA | 32 | 42 | 0.5067300 | MULTIPOLYGON (((10.42773 43… |
| PORDENONE | 29 | 40 | 0.4592241 | MULTIPOLYGON (((12.52112 46… |
| LATINA | 28 | 30 | 0.4433888 | MULTIPOLYGON (((12.94157 41… |
| CATANZARO | 26 | 28 | 0.4117181 | MULTIPOLYGON (((16.62624 39… |
| PESARO E URBINO | 26 | 27 | 0.4117181 | MULTIPOLYGON (((12.76836 43… |
| ROVIGO | 26 | 35 | 0.4117181 | MULTIPOLYGON (((12.32842 45… |
| FERRARA | 25 | 39 | 0.3958828 | MULTIPOLYGON (((12.01548 44… |
| MODENA | 25 | 27 | 0.3958828 | MULTIPOLYGON (((11.08931 44… |
| AREZZO | 24 | 33 | 0.3800475 | MULTIPOLYGON (((11.72191 43… |
| CREMONA | 24 | 26 | 0.3800475 | MULTIPOLYGON (((9.53991 45…. |
| ENNA | 24 | 63 | 0.3800475 | MULTIPOLYGON (((14.57105 37… |
| MANTOVA | 24 | 27 | 0.3800475 | MULTIPOLYGON (((10.66871 45… |
| ANCONA | 23 | 30 | 0.3642122 | MULTIPOLYGON (((13.50006 43… |
| ASCOLI PICENO | 21 | 25 | 0.3325416 | MULTIPOLYGON (((13.78665 43… |
| GORIZIA | 20 | 21 | 0.3167063 | MULTIPOLYGON (((13.53041 45… |
| VIBO VALENTIA | 20 | 22 | 0.3167063 | MULTIPOLYGON (((15.90041 38… |
| MONZA E DELLA BRIANZA | 19 | 28 | 0.3008709 | MULTIPOLYGON (((9.253944 45… |
| MASSA-CARRARA | 17 | 21 | 0.2692003 | MULTIPOLYGON (((9.903351 44… |
| PIACENZA | 16 | 21 | 0.2533650 | MULTIPOLYGON (((9.917413 45… |
| CASERTA | 15 | 16 | 0.2375297 | MULTIPOLYGON (((14.1601 41…. |
| TERNI | 15 | 15 | 0.2375297 | MULTIPOLYGON (((12.1291 42…. |
| BIELLA | 14 | 15 | 0.2216944 | MULTIPOLYGON (((8.129534 45… |
| NOVARA | 13 | 14 | 0.2058591 | MULTIPOLYGON (((8.496885 45… |
| PISTOIA | 13 | 14 | 0.2058591 | MULTIPOLYGON (((10.64794 44… |
| LODI | 11 | 12 | 0.1741884 | MULTIPOLYGON (((9.442017 45… |
| PRATO | 11 | 12 | 0.1741884 | MULTIPOLYGON (((11.1729 44…. |
| CAMPOBASSO | 8 | 8 | 0.1266825 | MULTIPOLYGON (((14.84675 42… |
| REGGIO NELL’EMILIA | 7 | 8 | 0.1108472 | MULTIPOLYGON (((10.72925 44… |
| BARLETTA-ANDRIA-TRANI | 5 | 5 | 0.0791766 | MULTIPOLYGON (((15.95552 41… |
| BENEVENTO | 4 | 4 | 0.0633413 | MULTIPOLYGON (((15.0303 41…. |
| CROTONE | 4 | 11 | 0.0633413 | MULTIPOLYGON (((17.05156 39… |
| FERMO | 4 | 4 | 0.0633413 | MULTIPOLYGON (((13.74288 43… |
| RIMINI | 4 | 4 | 0.0633413 | MULTIPOLYGON (((12.49431 44… |
9 Evolutionary Distinct and Globally Endangered
Όπως έχουμε αναφέρει σε προηγούμενο εργαστήριο, στο φυλογενετικό επίπεδο, οι ισοδύναμοι των παραδοσιακών, ταξινομικών δεικτών για τον εντοπισμό θερμών σημείων βιοποικιλότητας (ΘΣΒ), είναι η φυλογενετική ποικιλότητα [ΦΠ, sensu Faith (1992) – ισοδύναμη του αριθμού των ειδών], ο φυλογενετικός ενδημισμός [ΦΕ, sensu Rosauer et al. (2009) – ισοδύναμος του αριθμού των ενδημικών/ποσοστού ενδημισμού] και ο δείκτης EDGE [Evolutionary Distinct and Globally Endangered – EDGE, sensu Isaac et al. (2007) – ισοδύναμος της πιθανότητας εξαφάνισης]. Οι εν λόγω δείκτες ποσοτικοποιούν διαφορετικές πτυχές της εξελικτικής ποικιλότητας (Tucker et al., 2017).
Η ΦΠ είναι το άθροισμα του μήκους των κλάδων που συνδέουν ένα σύμπλεγμα ειδών με την ρίζα ενός φυλογενετικού δένδρου (Faith, 1992).
Ο ΦΕ ποσοτικοποιεί τον βαθμό στον οποίο ένα σημαντικό ποσοστό της ΦΠ εντοπίζεται αποκλειστικά εντός της εκάστοτε περιοχής μελέτης (Rosauer et al., 2009).
Ο δείκτης EDGE συνδυάζει την εξελικτική διακριτότητα (ED, η φυλογενετική απομόνωση ενός είδους) με το καθεστώς κινδύνου εξαφάνισης σε παγκόσμιο επίπεδο (GE) σύμφωνα με τα πρότυπα της IUCN (Isaac et al., 2007).
Οι δείκτες ΦΕ και EDGE αντιπροσωπεύουν σταθμισμένες παραλλαγές της ΦΠ σε γεωγραφική κλίμακα και καθεστώς απειλής, αντίστοιχα (Daru et al., 2019).
Συνεπώς ο δείκτης EDGE μπορεί να αποβεί ιδιαίτερα χρήσιμος κατά τον εντοπισμό και την προτεραιοποίηση των θερμών σημείων βιοποικιλότητας, ανεξαρτήτως κλίμακας (Daru et al., 2019).
Σήμερα λοιπόν, θα υπολογίσουμε τον δείκτη EDGE για τα φυτικά taxa τα οποία απαντώνται στην Ιταλία, σύμφωνα με τα δεδομένα τα οποία είναι διαθέσιμα στην βάση δεδομένων GBIF. θα πρέπει να τονιστεί ότι ο δείκτης EDGE πρέπει να υπολογίζεται για τα ενδημικά taxa της περιοχής μελέτης και όχι για το σύνολο των ειδών που απαντώνται σε αυτήν.
9.1 Υπολογισμός του δείκτη EDGE
Ας υπολογίσουμε τον δείκτη EDGE για τα είδη τα οποία απαντώνται στην Ιταλία.
## ============================================================================##
## Load the tree
## ============================================================================##
italy_tree <- readRDS("Italian phylogenetic tree.rds")
## ============================================================================##
## See which taxa are missing
## ============================================================================##
setdiff(italy_iucn_eval$Taxon, italy_tree$tip.label)## [1] "Cycadeoidea_etrusca"
## ============================================================================##
## Create the columns needed for the estimation of the EDGE index
## ============================================================================##
italy_iucn_eval_edge <- italy_iucn_eval %>% mutate(CRB = ifelse(Category_CriteriaB ==
"LC or NT", "NT", Category_CriteriaB)) %>% filter(Taxon != "Cycadeoidea_etrusca") %>%
dplyr::select(Taxon, CRB) %>% rename(Redlists = CRB, Species = Taxon) %>% as.data.frame() %>%
dplyr::select(Species, Redlists)
## ============================================================================##Εάν σας δείτε κάποιο error στην επόμενη εντολή, τρέξτε τις σειρές του κώδικα που εμφανίζονται ως σχόλια.
## ============================================================================##
## Estimate EDGE -----
## ============================================================================##
italian_edge <- phyloregion::EDGE(italy_iucn_eval_edge, italy_tree, Redlist = "Redlists",
species = "Species")
## ============================================================================##
## ============================================================================##
## If the above throws an error, run the following chain of commands
## ============================================================================##
## italian_edge <- phyloregion::evol_distinct(italy_tree, type =
## 'fair.proportion') %>% enframe() %>% dplyr::rename(Taxon = name, ED = value)
## %>% inner_join(., italy_iucn_eval) %>% mutate(CRB = ifelse(Category_CriteriaB
## == 'CR', 0.999, ifelse(Category_CriteriaB == 'EN', 0.67,
## ifelse(Category_CriteriaB == 'VU', 0.1, 0.01))), EDGE = EDGE_fnc(ED, CRB)) %>%
## dplyr::select(Taxon, ED, EDGE, everything()) map_edge <- italian_edge$EDGE
## names(map_edge) <- italian_edge$Taxon
## ============================================================================#### ============================================================================##
## Create a duplicate for later analysis
## ============================================================================##
italian_edge_tbl <- italian_edge %>% enframe() %>% dplyr::rename(Taxon = name, EDGE = value) %>%
inner_join(., clean_dts) %>% dplyr::select(order, family, genus, Taxon, EDGE)
## ============================================================================##Ας δούμε ποια είναι η κατανομή των τιμών του δείκτη EDGE για τα φυτικά taxa της Ιταλίας, κάποια στατιστικά στοιχεία για τον δείκτη αυτόν, καθώς και ποια είδη (αλλά και γένη και οικογένειες) εμφανίζουν τις υψηλότερες τιμές για τον δείκτη EDGE και κατά συνέπεια θα πρέπει να δώσουμε μεγαλύτερη έμφαση στην προτεραιοποίηση τους:
## ============================================================================##
## Histogram
## ============================================================================##
hist(italian_edge_tbl$EDGE, main = "Distribution of EDGE values in Italy", xlab = "EDGE")## ============================================================================##
## Basic statistical features
## ============================================================================##
skimr::skim(italian_edge_tbl)| Name | italian_edge_tbl |
| Number of rows | 6314 |
| Number of columns | 5 |
| _______________________ | |
| Column type frequency: | |
| character | 4 |
| numeric | 1 |
| ________________________ | |
| Group variables | None |
Variable type: character
| skim_variable | n_missing | complete_rate | min | max | empty | n_unique | whitespace |
|---|---|---|---|---|---|---|---|
| order | 0 | 1 | 6 | 15 | 0 | 59 | 0 |
| family | 0 | 1 | 7 | 16 | 0 | 179 | 0 |
| genus | 0 | 1 | 3 | 20 | 0 | 1222 | 0 |
| Taxon | 0 | 1 | 8 | 48 | 0 | 6314 | 0 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| EDGE | 0 | 1 | 3.02 | 0.86 | 0.64 | 2.4 | 2.94 | 3.57 | 6.32 | ▁▇▇▂▁ |
## ============================================================================##
## Mean EDGE values by genus
## ============================================================================##
italian_edge_tbl %>% group_by(family, genus) %>% summarise(Mean_EDGE = mean(EDGE),
Species_Number = sum(NROW(genus))) %>% arrange(desc(Species_Number)) %>% ungroup() %>%
head()| family | genus | Mean_EDGE | Species_Number |
|---|---|---|---|
| Orchidaceae | Ophrys | 3.052281 | 122 |
| Cyperaceae | Carex | 2.166766 | 97 |
| Asteraceae | Centaurea | 2.121417 | 93 |
| Asteraceae | Hieracium | 2.114490 | 84 |
| Euphorbiaceae | Euphorbia | 2.933861 | 70 |
| Caryophyllaceae | Silene | 2.988936 | 66 |
## ============================================================================##
## Find out which species have the top-1% of EDGE values
## ============================================================================##
italian_edge_tbl %>% arrange(desc(EDGE)) %>% dplyr::select(Taxon, EDGE) %>% filter(EDGE >=
quantile(italian_edge_tbl$EDGE, probs = 0.99, na.rm = T)) %>% head()| Taxon | EDGE |
|---|---|
| Huperzia_selago | 6.315989 |
| Spinulum_annotinum | 6.087944 |
| Equisetum_hyemale_subsp_hyemale | 6.055210 |
| Equisetum_meridionale | 6.055210 |
| Cycas_revoluta | 6.036176 |
| Equisetum_arvense | 5.789793 |
9.2 Χωρικά πρότυπα του δείκτη EDGE - Θερμά σημεία βιοποικιλότητας
Τώρα μπορούμε να υπολογίσουμε ποιες περιοχές της Ιταλίας εμφανίζουν υψηλές τιμές για τον δείκτη EDGE και μπορούν να θεωρηθούν ως θερμά ή ψυχρά σημεία βιοποικιλότητας. Πρώτα θα χρειαστεί να εισάγαγουμε το αρχείο το οποίο είχαμε δημιουργήσει σε προηγούμενο εργαστήριο και αφορά τη μήτρα παρουσίας/απουσίας των φυτικών taxa στην Ιταλία:
## ============================================================================##
## Load the sparse matrix
## ============================================================================##
pt <- readRDS("Sparse matrix and relevant polygons for Italy.rds")
## ============================================================================##Και στη συνέχεια, μπορούμε να δημιουργήσουμε ένα χωρικό αντικείμενο το οποίο θα περιέχει πληροφορία σχετικά με το ποιες περιοχές δρουν ως θερμά ή ψυχρά σημεία βιοποικιλότητας (ταξινομικής και φυλογενετικής), καθώς και ως θερμά ή ψυχρά σημεία κινδυνού εξαφάνισης (δείκτης EDGE):
##============================================================================##
## Create a sf object with hot-/cold-spots info
##============================================================================##
hot_cold_spots <- pt$poly_shp %>%
## Convert to sf
st_as_sf() %>%
## Create the hot- and cold-spots
mutate(Cold_SR = coldspots(richness,
prob = 5),
Hot_SR = hotspots(richness,
prob = 5),
## Corrected weighted endemism
Cold_CWE = coldspots(weighted_endemism(pt$comm_dat),
prob = 5),
Hot_CWE = hotspots(weighted_endemism(pt$comm_dat),
prob = 5),
## Phylogenetic diversity
Cold_PD = coldspots(PD(match_phylo_comm(italy_tree, pt$comm_dat)$com,
match_phylo_comm(italy_tree, pt$comm_dat)$phy),
prob = 5),
Hot_PD = hotspots(PD(match_phylo_comm(italy_tree, pt$comm_dat)$com,
match_phylo_comm(italy_tree, pt$comm_dat)$phy),
prob = 5),
## Phylogenetic endemism
Cold_PE = coldspots(phylo_endemism(match_phylo_comm(italy_tree, pt$comm_dat)$com,
match_phylo_comm(italy_tree, pt$comm_dat)$phy),
prob = 5),
Hot_PE = hotspots(phylo_endemism(match_phylo_comm(italy_tree, pt$comm_dat)$com,
match_phylo_comm(italy_tree, pt$comm_dat)$phy),
prob = 5),
## EDGE
EDGE = map_trait(pt$comm_dat,
## If an error occured previously, change that with map_edge
italian_edge,
FUN = median,
shp = pt$poly_shp)$traits,
Cold_EDGE = coldspots(EDGE, prob = 5),
Hot_EDGE = hotspots(EDGE, prob = 5))9.3 Οπτικοποίηση
Ας δημιουργήσουμε ένα νέο αρχείο το οποίο θα χρησιμοποιήσουμε για την οπτικοποίηση των αποτελέσματων μας και το οποίο θα περιέχει μια νέα μεταβλητή, η οποία θα υποδεικνύει ποια κελιά αποτελούν θερμά ή ψυχρά σημεία βιοποικιλότητας για όλους τους δείκτες που χρησιμοποιήσαμε:
## ============================================================================##
## Create a new variable
## ============================================================================##
sp_hot_cold_spots <- hot_cold_spots %>% as_Spatial()
## ============================================================================##
## ============================================================================##
## Quick plot
## ============================================================================##
plot(sp_hot_cold_spots, border = "grey", col = "lightgrey", main = "Diversity Hotspots and Coldspots")
plot(sp_hot_cold_spots[(sp_hot_cold_spots@data$Hot_EDGE == 1), ], col = "red", add = TRUE,
border = NA)
plot(sp_hot_cold_spots[(sp_hot_cold_spots@data$Cold_EDGE == 1), ], col = "blue",
add = TRUE, border = NA)
legend("bottomleft", fill = c("red", "blue"), legend = c("hotspots", "coldspots"),
bty = "n", inset = 0.092)10 Appendix: R code
##===========================================================================##
## Do NOT run this code
##===========================================================================##
## Global options
options(Encoding="UTF-8")
library(knitr)
library(rmdformats)
## Global options
options(max.print="75")
opts_chunk$set(echo=TRUE,
cache=TRUE,
prompt=FALSE,
tidy=TRUE,
# comment=NA,
message=FALSE,
warning=FALSE,
eval = TRUE)
opts_knit$set(width=75)
##===========================================================================##
htmltools::img(src = knitr::image_uri(file.path("G:/Academic/Lessons/EMB/Labs/Labs/TRUE RMD Files/HBS Seminar/Logo_HBS/ebe-LOGO-005-300_1_70.jpg")),
alt = 'logo',
style = 'position:absolute; top:0; right:0; padding:5px;')
##===========================================================================##
## Load them
##===========================================================================##
library(raster)
library(sf)
library(tidyverse)
library(ConR)
library(red)
library(rasterVis)
library(ggspatial)
library(picante)
library(caper)
library(phyloregion)
library(ggsflabel)
library(cowplot)
##===========================================================================##
##===========================================================================##
## Load the spatial data for Italian provinces
##===========================================================================##
Italy <- read_sf('Italin provinces.shp') %>%
st_transform(4326) %>%
as_Spatial()
##===========================================================================##
##===========================================================================##
## Load the spatial data for Italian counties
##===========================================================================##
Italian_counties <- read_sf("reg2011_g.shp") %>%
mutate(NOME_REG = str_replace_all(NOME_REG,
c("VALLE D'AOSTA/VALLÉE D'AOSTE\r\nVALLE D'AOSTA/VALLÉE D'AOSTE" = "AOSTA",
"TRENTINO-ALTO ADIGE/SUDTIROL" = "TRENTINO",
"FRIULI VENEZIA GIULIA" ="FRIULI",
"EMILIA-ROMAGNA" ="EMILIA"))) %>%
st_transform(4326)
##===========================================================================##
##===========================================================================##
## Keep only the names in Trentino
##===========================================================================##
Italy_sf_tr <- Italy %>%
st_as_sf() %>%
st_join(Italian_counties) %>%
filter(str_detect(NOME_REG, 'TREN'))
##===========================================================================##
##============================================================================##
## Load the biotic cleaned data
##============================================================================##
data_it_cleaned <- readRDS('Italian data cleaned.rds')
##============================================================================##
##============================================================================##
## First load the excel file with the correct species names
##============================================================================##
nms <- readxl::read_excel('italian_names.xlsx')
##============================================================================##
##============================================================================##
## Keep only one instance per taxon from our occurrence data
##============================================================================##
nms_flags <- data_it_cleaned %>%
distinct(scientificName) %>%
mutate(Taxon = nms$scientificName,
action = nms$action) %>%
as_tibble()
##============================================================================##
##============================================================================##
## Then create the coords for the correct taxa
##============================================================================##
coords <- data_it_cleaned %>%
full_join(nms_flags) %>%
as_tibble() %>%
filter(!str_detect(action, 'DEL')) %>%
dplyr::select(decimalLongitude,
decimalLatitude,
Taxon)
clean_dts <- data_it_cleaned %>%
full_join(nms_flags) %>%
as_tibble() %>%
filter(!str_detect(action, 'DEL')) %>%
distinct(Taxon, .keep_all = T)
##============================================================================##
pre {
max-height: 400px;
overflow-y: auto;
}
pre[class] {
max-height: 400px;
}
library(magick)
nms <- list.files('MAGICK',
full.names = T)
lapply(nms, image_read) %>%
image_join() %>%
image_animate(fps = 1)
aegean %>%
download_this(
output_name = "PVA",
output_extension = ".xlsx",
button_label = "Download the dataset",
button_type = "success"
)
library(downloadthis)
aegean <- readxl::read_excel("./PVA_Crete.xlsx")
##==================================================================================##
## Load the data
##==================================================================================##
cretan_taxa <- readxl::read_excel("./PVA_Crete.xlsx") %>%
dplyr::rename(decimalLongitude = lon,
decimalLatitude = lat)
cretan_taxa$Taxon <- factor(cretan_taxa$Taxon)
levels(cretan_taxa$Taxon)
##==================================================================================##
## Create 3 objects
##==================================================================================##
all_dila <- cretan_taxa %>% filter(Taxon == 'Allium dilatatum')
camp_jac <- cretan_taxa %>% filter(Taxon == 'Campanula jacquinii subsp. jacquinii')
tul_doer <- cretan_taxa %>% filter(Taxon == 'Tulipa doerfleri')
##==================================================================================##
##==================================================================================##
## Convert to matrices
##==================================================================================##
all_dila_mat <- all_dila %>% dplyr::select(decimalLongitude:decimalLatitude) %>% as.matrix()
camp_jac_mat <- camp_jac %>% dplyr::select(decimalLongitude:decimalLatitude) %>% as.matrix()
tul_doer_mat <- tul_doer %>% dplyr::select(decimalLongitude:decimalLatitude) %>% as.matrix()
##==================================================================================##
##==================================================================================##
## Area of occupancy
##==================================================================================##
aoo(all_dila_mat)
aoo(camp_jac_mat)
aoo(tul_doer_mat)
##==================================================================================##
##==================================================================================##
## Extent of Occurrence
##==================================================================================##
eoo(all_dila_mat)
eoo(camp_jac_mat)
eoo(tul_doer_mat)
##==================================================================================##
##==================================================================================##
## Load the rds files
##==================================================================================##
allium_future <- readRDS('./RDS/Allium future.rds')
campanula_future <- readRDS('./RDS/Campanula future.rds')
tulipa_future <- readRDS('./RDS/Tulipa future.rds')
##==================================================================================##
##==================================================================================##
## Area of occupancy
##==================================================================================##
aoo(allium_future)
aoo(campanula_future)
aoo(tulipa_future)
##==================================================================================##
##==================================================================================##
## Extent of Occurrence
##==================================================================================##
eoo(allium_future)
eoo(campanula_future)
eoo(tulipa_future)
##==================================================================================##
##============================================================================##
## First, let's select only the required columns
##============================================================================##
data_it_sel <- coords %>%
dplyr::select(decimalLatitude,
decimalLongitude,
Taxon)
##============================================================================##
##============================================================================##
## Then, let's evaluate our taxa
##============================================================================##
italy_iucn_eval <- IUCN.eval(data_it_sel,
## Here we use the spatial object for the country
## we want
country_map = Italy,
## We can run the function in parallel to save
## time
parallel = T,
## Change the number of cores accordingly
NbeCores = 23 ) %>%
## Here we create a new column (EOO_tr) and we substitute the NA values with
## the values from the AOO column (Why?)
mutate(EOO_tr = ifelse(is.na(EOO),
AOO,
EOO)) %>%
as_tibble() %>%
dplyr::select(taxa,
EOO_tr,
everything())
##============================================================================##
## Let's see the result
##============================================================================##
italy_iucn_eval
##============================================================================##
## First we run the convex-hull method
##============================================================================##
eooshp_convexhull <- EOO.computing(data_gr_sel,
exclude.area = T,
country_map = Greece,
export_shp = T,
write_shp = T,
Name_Sp = 'scientificName',
parallel = T,
NbeCores = 4, ## Change them
show_progress = T,
method.less.than3 = 'arbitrary') ## What does this option mean?
convexhull <- eooshp_convexhull$spatial.polygon_1
##============================================================================##
##============================================================================##
## Then the alpha-hull method
##============================================================================##
eooshp_ahull <- EOO.computing(data_gr_sel,
exclude.area = T,
country_map = Greece,
export_shp = T,
write_shp = T,
Name_Sp = 'scientificName',
parallel = T,
NbeCores = 4, ## Change them
show_progress = T,
method.less.than3 = 'arbitrary',
method.range = 'alpha.hull')
ahull <- eooshp_ahull$spatial.polygon_1
##============================================================================##
##============================================================================##
## First create a common column and then join the occurrences with the
## evaluations
##============================================================================##
italy_iucn_eval <- italy_iucn_eval %>% rename(Taxon = taxa)
data_joined <- full_join(coords, italy_iucn_eval)
##============================================================================##
##============================================================================##
## Then select only the endangered taxa (CR, EN, VU)
##============================================================================##
endangered <- data_joined %>% filter(!str_detect(Category_code, 'LC'))
##============================================================================##
##============================================================================##
## Convert Italy to an sf spatial object
##============================================================================##
Italy_sf <- Italy %>% st_as_sf()
##============================================================================##
##============================================================================##
## Convert the endangered species to a sf spatial object
##============================================================================##
data_it_sf <- st_as_sf(endangered,
coords = c("decimalLongitude",
"decimalLatitude"),
crs = 4326)
##============================================================================##
##============================================================================##
## Make sure that all these occurrences lie within Greece
##============================================================================##
point_endangered <- data_it_sf[Italy_sf,]
##============================================================================##
##============================================================================##
## Intersect Greece with these occurrences
##============================================================================##
endangered_intersection <- Italy_sf %>%
st_join(point_endangered) %>%
group_by(NOME_PRO) %>% ## Here we group the data based on the county names
## First we create a variable containing the sum of the species in each county
summarize(Species_Number = sum(NROW(unique(Taxon))),
## Then we just sum the number of occurrences in each county
Occurrence_Number = n(),
## Finally, we calculate the proportion of the endangered species in
## each county
Risk = sum(NROW(unique(Taxon)))/sum(NROW(unique(data_joined$Taxon)))*100
)
##============================================================================##
##============================================================================##
## Colour function
##============================================================================##
col2alpha <- function(col, alpha) {
col_rgb <- col2rgb(col)/255
rgb(col_rgb[1], col_rgb[2], col_rgb[3],
alpha = alpha)
}
##============================================================================##
##============================================================================##
## Create the main map
##============================================================================##
main_map <- ggplot(Italy_sf) +
geom_sf(fill = "antiquewhite1") +
geom_sf() +
geom_sf(data = Italy_sf,
fill = NA) +
geom_sf(aes(fill = Risk),
data = endangered_intersection,
inherit.aes = FALSE) +
theme(axis.title = element_blank(),
legend.position = "bottom") +
coord_sf(expand = FALSE) +
scale_fill_viridis_c(option = "viridis",
trans = "sqrt",
limits = c(0, 12.4),
breaks = c(0.4, 3, 12.3),
labels = scales::label_percent(scale = 1),
name = expression(bold('Extinction risk'))) +
geom_rect(
xmin = st_bbox(Italy_sf_tr)[[1]],
ymin = st_bbox(Italy_sf_tr)[[2]],
xmax = st_bbox(Italy_sf_tr)[[3]],
ymax = st_bbox(Italy_sf_tr)[[4]],
fill = NA,
colour = "red",
size = 0.6
) +
annotation_scale(location = "bl",
width_hint = 0.25) +
annotation_north_arrow(location = "bl",
which_north = "true",
pad_x = unit(0.75, "in"),
pad_y = unit(0.5, "in"),
style = north_arrow_fancy_orienteering) +
geom_sf(data = Italian_counties,
fill = NA,
lwd = 1,
color = 'black') +
geom_sf_label_repel(aes(label = NOME_REG),
data = Italian_counties) +
theme(panel.grid.major = element_line(color = gray(.5),
linetype = "dashed",
size = 0.5),
panel.background = element_rect(fill = col2alpha('steelblue', 0.2))) +
labs(title = 'Extinction risk in Italian counties',
subtitle = 'Source: GBIF')
##============================================================================##
##============================================================================##
## Create the inset map
##============================================================================##
side_map <- ggplot(Italy_sf) +
geom_sf(fill = "antiquewhite1") +
geom_sf() +
geom_sf(data = Italy_sf,
fill = NA) +
geom_sf(aes(fill = Risk),
data = endangered_intersection,
inherit.aes = FALSE) +
theme(axis.title = element_blank(),
legend.position = "bottom") +
coord_sf(expand = FALSE) +
scale_fill_viridis_c(option = "viridis",
trans = "sqrt",
limits = c(0, 12.4),
breaks = c(0.4, 3, 12.3),
labels = scales::label_percent(scale = 1),
name = 'Extinction risk') +
geom_rect(
xmin = st_bbox(Italy_sf_tr)[[1]],
ymin = st_bbox(Italy_sf_tr)[[2]],
xmax = st_bbox(Italy_sf_tr)[[3]],
ymax = st_bbox(Italy_sf_tr)[[4]],
fill = NA,
colour = "red",
size = 0.6
) +
geom_sf(data = Italian_counties,
fill = NA,
lwd = 1,
color = 'black') +
geom_sf_label_repel(aes(label = NOME_PRO),
data = Italy_sf_tr) +
coord_sf(xlim = c(st_bbox(Italy_sf_tr)[[1]],
st_bbox(Italy_sf_tr)[[3]]),
ylim = c(st_bbox(Italy_sf_tr)[[2]],
st_bbox(Italy_sf_tr)[[4]]),
expand = FALSE) +
theme_void() +
theme(legend.position = "none")
##============================================================================##
##============================================================================##
## And the final map
##============================================================================##
main_map %>%
ggdraw() +
draw_plot(side_map,
# The distance along a (0,1) x-axis to draw the left edge of the plot
x = 0.625,
# The distance along a (0,1) y-axis to draw the bottom edge of the plot
y = 0.65,
# The width and height of the plot expressed as proportion of the entire
# ggdraw object
width = 0.35,
height = 0.35)
##============================================================================##
## The solution
##============================================================================##
endangered_intersection %>% arrange(desc(Risk))
##============================================================================##
## Load the tree
##============================================================================##
italy_tree <- readRDS('Italian phylogenetic tree.rds')
##============================================================================##
## See which taxa are missing
##============================================================================##
setdiff(italy_iucn_eval$Taxon, italy_tree$tip.label)
##============================================================================##
## Create the columns needed for the estimation of the EDGE index
##============================================================================##
italy_iucn_eval_edge <- italy_iucn_eval %>%
mutate(CRB = ifelse(Category_CriteriaB == 'LC or NT',
'NT',
Category_CriteriaB)) %>%
filter(Taxon != 'Cycadeoidea_etrusca') %>%
dplyr::select(Taxon, CRB) %>%
rename(Redlists = CRB,
Species = Taxon) %>%
as.data.frame() %>%
dplyr::select(Species, Redlists)
##============================================================================##
italian_edge <- readRDS('Italian EDGE.rds')
##============================================================================##
## Estimate EDGE -----
##============================================================================##
italian_edge <- phyloregion::EDGE(italy_iucn_eval_edge,
italy_tree,
Redlist = "Redlists",
species = "Species")
##============================================================================##
##============================================================================##
## If the above throws an error, run the following chain of commands
##============================================================================##
# italian_edge <- phyloregion::evol_distinct(italy_tree,
# type = 'fair.proportion') %>%
# enframe() %>%
# dplyr::rename(Taxon = name,
# ED = value) %>%
# inner_join(., italy_iucn_eval) %>%
# mutate(CRB = ifelse(Category_CriteriaB == 'CR', 0.999,
# ifelse(Category_CriteriaB == 'EN', 0.67,
# ifelse(Category_CriteriaB == 'VU', 0.1, 0.01))),
# EDGE = EDGE_fnc(ED, CRB)) %>%
# dplyr::select(Taxon, ED, EDGE, everything())
#
# map_edge <- italian_edge$EDGE
# names(map_edge) <- italian_edge$Taxon
##============================================================================##
##============================================================================##
## Create a duplicate for later analysis
##============================================================================##
italian_edge_tbl <- italian_edge %>%
enframe() %>%
dplyr::rename(Taxon = name,
EDGE = value) %>%
inner_join(., clean_dts) %>%
dplyr::select(order, family, genus, Taxon, EDGE)
##============================================================================##
##============================================================================##
## Histogram
##============================================================================##
hist(italian_edge_tbl$EDGE,
main = 'Distribution of EDGE values in Italy',
xlab = 'EDGE')
##============================================================================##
## Basic statistical features
##============================================================================##
skimr::skim(italian_edge_tbl)
##============================================================================##
## Mean EDGE values by genus
##============================================================================##
italian_edge_tbl %>%
group_by(family, genus) %>%
summarise(Mean_EDGE = mean(EDGE),
Species_Number = sum(NROW(genus))) %>%
arrange(desc(Species_Number)) %>%
ungroup() %>%
head()
##============================================================================##
## Find out which species have the top-1% of EDGE values
##============================================================================##
italian_edge_tbl %>%
arrange(desc(EDGE)) %>%
dplyr::select(Taxon, EDGE) %>%
filter(EDGE >= quantile(italian_edge_tbl$EDGE,
probs = 0.99,
na.rm = T)) %>%
head()
##============================================================================##
## Load the sparse matrix
##============================================================================##
pt <- readRDS('Sparse matrix and relevant polygons for Italy.rds')
##============================================================================##
##============================================================================##
## Create a sf object with hot-/cold-spots info
##============================================================================##
hot_cold_spots <- pt$poly_shp %>%
## Convert to sf
st_as_sf() %>%
## Create the hot- and cold-spots
mutate(Cold_SR = coldspots(richness,
prob = 5),
Hot_SR = hotspots(richness,
prob = 5),
## Corrected weighted endemism
Cold_CWE = coldspots(weighted_endemism(pt$comm_dat),
prob = 5),
Hot_CWE = hotspots(weighted_endemism(pt$comm_dat),
prob = 5),
## Phylogenetic diversity
Cold_PD = coldspots(PD(match_phylo_comm(italy_tree, pt$comm_dat)$com,
match_phylo_comm(italy_tree, pt$comm_dat)$phy),
prob = 5),
Hot_PD = hotspots(PD(match_phylo_comm(italy_tree, pt$comm_dat)$com,
match_phylo_comm(italy_tree, pt$comm_dat)$phy),
prob = 5),
## Phylogenetic endemism
Cold_PE = coldspots(phylo_endemism(match_phylo_comm(italy_tree, pt$comm_dat)$com,
match_phylo_comm(italy_tree, pt$comm_dat)$phy),
prob = 5),
Hot_PE = hotspots(phylo_endemism(match_phylo_comm(italy_tree, pt$comm_dat)$com,
match_phylo_comm(italy_tree, pt$comm_dat)$phy),
prob = 5),
## EDGE
EDGE = map_trait(pt$comm_dat,
## If an error occured previously, change that with map_edge
italian_edge,
FUN = median,
shp = pt$poly_shp)$traits,
Cold_EDGE = coldspots(EDGE, prob = 5),
Hot_EDGE = hotspots(EDGE, prob = 5))
##============================================================================##
## Create a new variable
##============================================================================##
sp_hot_cold_spots <- hot_cold_spots %>%
as_Spatial()
##============================================================================##
##============================================================================##
## Quick plot
##============================================================================##
plot(sp_hot_cold_spots,
border = "grey",
col = "lightgrey",
main = "Diversity Hotspots and Coldspots")
plot(sp_hot_cold_spots[(sp_hot_cold_spots@data$Hot_EDGE == 1), ],
col = "red",
add = TRUE,
border = NA)
plot(sp_hot_cold_spots[(sp_hot_cold_spots@data$Cold_EDGE == 1), ],
col = "blue",
add = TRUE,
border = NA)
legend("bottomleft",
fill = c("red", "blue"),
legend = c("hotspots", "coldspots"),
bty = "n",
inset = .092)
##===========================================================================##
## Do NOT run this code
##===========================================================================##
labs = knitr::all_labels(echo = TRUE)
##===========================================================================##Βιβλιογραφία
Bachman, S.P., Nic Lughadha, E.M., & Rivers, M.C. (2018) Quantifying progress toward a conservation assessment for all plants. Conservation biology, 32, 516–524.
Daru, B.H., Roux, P.C. le, Gopalraj, J., Park, D.S., Holt, B.G., & Greve, M. (2019) Spatial overlaps between the global protected areas network and terrestrial hotspots of evolutionary diversity. Global Ecology and Biogeography, 28, 757–766.
Dauby, G., Stévart, T., Droissart, V., Cosiaux, A., Deblauwe, V., Simo-Droissart, M., Sosef, M.S., Lowry, P.P., Schatz, G.E., Gereau, R.E., & others (2017) ConR: An r package to assist large-scale multispecies preliminary conservation assessments using distribution data. Ecology and evolution, 7, 11292–11303.
Faith, D.P. (1992) Conservation evaluation and phylogenetic diversity. Biological conservation, 61, 1–10.
Isaac, N.J., Turvey, S.T., Collen, B., Waterman, C., & Baillie, J.E. (2007) Mammals on the edge: Conservation priorities based on threat and phylogeny. PloS one, 2, e296.
Kougioumoutzis, K., Kokkoris, I.P., Panitsa, M., Trigas, P., Strid, A., & Dimopoulos, P. (2020a) Plant diversity patterns and conservation implications under climate-change scenarios in the mediterranean: The case of crete (aegean, greece). Diversity, 12, 270.
Kougioumoutzis, K., Kokkoris, I.P., Panitsa, M., Trigas, P., Strid, A., & Dimopoulos, P. (2020b) Spatial phylogenetics, biogeographical patterns and conservation implications of the endemic flora of crete (aegean, greece) under climate change scenarios. Biology, 9, 199.
Rosauer, D., Laffan, S.W., Crisp, M.D., Donnellan, S.C., & Cook, L.G. (2009) Phylogenetic endemism: A new approach for identifying geographical concentrations of evolutionary history. Molecular ecology, 18, 4061–4072.
Stévart, T., Dauby, G., Lowry, P., Blach-Overgaard, A., Droissart, V., Harris, D., Mackinder, B., Schatz, G., Sonké, B., Sosef, M., & others (2019) A third of the tropical african flora is potentially threatened with extinction. Science Advances, 5, eaax9444.
Tucker, C.M., Cadotte, M.W., Carvalho, S.B., Davies, T.J., Ferrier, S., Fritz, S.A., Grenyer, R., Helmus, M.R., Jin, L.S., Mooers, A.O., & others (2017) A guide to phylogenetic metrics for conservation, community ecology and macroecology. Biological Reviews, 92, 698–715.